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ABSTRACT 

We present ALMA (Cycle 0) band-6 and band-3 observations of the transition disk Sz 91. The 
disk inclination and position angle are determined to be i = 49.5°±3.5° and PA = 18.2°±3.5° and 
the dusty and gaseous disk are detected up to ^ 220 au and ^ 400 au from the star, respectively. 
Most importantly, our continuum observations indicate that the cavity size in the mm-sized dust 
distribution must be ^ 97 au in radius, the largest cavity observed around a T Tauri star. Our 
data clearly confirms the presence of (2-1) well inside the dust cavity. Based on these 

observational constrains we developed a disk model that simultaneously accounts for the 
and continuum observations (i.e., gaseous and dusty disk). According to our model, most of the 
millimeter emission comes from a ring located between 97 and 140 au. We also find that the 
dust cavity is divided into an innermost region largely depleted of dust particles ranging from the 
dust sublimation radius up to 85 au, and a second, moderately dust-depleted region, extending 
from 85 to 97 au. The extremely large size of the dust cavity, the presence of gas and small dust 
particles within the cavity and the accretion rate of Sz 91 are consistent with the formation of 
multiple (giant) planets. 

Subject headings: planet-disk interactions—protoplanetary disks—stars: variables: T Tauri, Herbig 
Ae/Be —stars: individual (Sz 91) 
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1. Introduction 


Simulations of disk-planet interactions show 
that a giant planet embedded in a circumstellar 
disk should open a gap or carve an inner hole in 
the disk ( Lin fc Papaloizou~||1986 1. These gaps or 
holes are imprinted in the spectral energy distri¬ 
butions (SEDs) of young stellar objects as reduced 
near and/or mid infrared excess emission. Proto¬ 
planetary disks showing this observational feature 
are usually called “transition disks”, although dif¬ 
ferent names and classihcations are found in the 
literature ( Evans et al. [2009 ). Apart from planet 
formation, three alternative mechanisms can cre¬ 
ate opacity holes in the inner disk: grain growth 
( Dullemond fc Dominik |[20Q5 ), photoevaporation 


(Alexander et al. 2006), and tidal truncation by 


close stellar companions (Artymowicz & Lubow 
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The latest models of photoevaporation show 
that only transition disks with small accretion 
rates and/or small inner cavities can be explained 


by this process (Owen et al. 20121. While grain 


growth and transport effects can explain the dips 
in the SEDs of transition disks, the most recent 
models fail in reproducing the large cavities ob¬ 
served in millimeter images of accreting transition 
disks (Birnstiel et al. ||2013|. Current planet for¬ 


mation models predict that a forming planet may 
reduce (or remove if the planet is massive enough) 
the amount of small particles within the disk’s cav¬ 
ity while allowing the gas to flow through the cav¬ 
ity (e.g., |Rice et al. ~ 2006[ Lubow & D’Angelo 


2006[ [Dodson-Robinson fc Salyk 1120111^0 et al. 


The transition disk Sz91, at 200 pc (Merin et 
~||2008 ), has been classified as a planet forming 


al. 


candidate based on its lack of detected compan¬ 
ions down to 30 au, its disk mass, accretion rate 
(7.9 X 1O“^^M0 yr“^), and SED shape (Romero et 


al. 2012). Recently, Tsukagoshi et al. (2014) 


(hereafter T2014) resolved the inner cavity in K- 
band polarized light with Hi CIAO/Subaru, and 
find evidences for a gap in the continuum at 345 
GHz (870 fj,m) from SMA observations. These au¬ 
thors also detect the outer gaseous disk using the 
^^CO(3—2) line. We here use ALMA cycle 0 band- 
6 and band-3 observations to describe the disk in 
more detail. Our observations reveal a huge cavity 
of ^97 au in radius in the continuum at 230 GHz 
(1.3 mm), and confirm the presence of ^^CO(2— 1) 
inside the cavity. We also derive stronger con¬ 
strains on the disk’s geometry and orientation. 


2. Observations and Data reduction 
2.1. Interferometry 

Sz91 was observed with ALMA band-6 (231 
GHz) and band-3 (110 GHz) during Gycle-0 (pro¬ 
gram 2011.0.00733, PI: M.R. Schreiber). The ob¬ 
servations are summarized in Table [TJ 

Sz 91 was observed in five different epochs with 
different weather conditions and antenna config¬ 
urations for the band-6 data. The system tem¬ 
perature ranged from 53 to 96 K. For the band- 
3 data, the ALMA correlator was configured to 
provide two spectral windows centered on the 
continuum, and two spectral windows centered 


on the ^^CO(l — 0) (110.20135 GHz) and the 
Gi®O(l- 0) (109.78218 GHz) lines. The total 
bandwidth of the band-3 observations was 7.5 
GHz (4X 1.875 GHz). For the band-6 observations, 
the ALMA correlator was initially configured to 
provide 1 spectral window centered on the con¬ 
tinuum, and three spectral windows centered on 
the ^^CO(2 — 1), 230.538 GHz (hereafter ), 
^^GO(2—1) and C^®0(2—1) lines. Unfortunately, 
because of technical problems only one sideband 
of the correlator could be configured, and only 
the spectral windows centered on the continuum 
and on the ^^GO line could be produced. The 
total bandwidth of the band-6 observations was 
2 X 1.875 GHz. Both ALMA bands were sampled 
at 0.488 MHz (e.g., 0.635 km/s in the ^^GO line). 

The quasar QSO B1730-130 and the AGN 
ICRF J160431.0-444131 were used as bandpass 
and primary phase calibrator in both bands, re¬ 
spectively. Neptune and Titan were used to cali¬ 
brate in amplitude in band-3 and band-6, respec¬ 
tively. The observations of the calibrators were 
alternated with the science observations. The in¬ 
dividual exposure times for the calibrators and the 
science targets was 6.05 seconds, amounting to a 
total exposure time for the science observations 
of 226 seconds (3.77 minutes) in band-6 and 24 
seconds (0.4 minutes) in band-3. Phase correction 
based on WVR measurements was performed in 
offline mode, as part of the basic ALMA Gycle 0 
corrections. We used the dispersion of the band- 
6 flux calibrations to estimate a flux uncertainty 
of « 15%. This value is a lower limit since it 
does not include uncertainties from the amplitude 
calibratiorH 


The observations were processed using the 
Gommon Astronomy Software Application (CASA) 
package (McMullin et al. 2007). The visibili¬ 
ties were Fourier transformed and deconvolved, 
using natural weighting, with the CLEAN algo¬ 
rithm (Hogbom, J. A. 1974). For the band-6 


data, we combined the four datasets to increase 
the uv coverage before cleaning. The cleaned im¬ 
age produced from this concatenated dataset had 
lower rms than the images created from individu¬ 
ally cleaned datasets. After combining the band-6 


^According to the documentation found in 
http://www.almaobservatory.org/ the amplitude cali¬ 
bration uncertainties are < 5%. 
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Table 1: Observing log. 


LO 

[GHz] 

Date 

[dd/mm/yy] 

Flux 

Galibrator 

Tv 

Zenith 

PWV 

[mm] 

texp 

[min] 

Tsys 

[K] 

bandwidth 

[GHz] 

231.46 

18/06/2012 

Titan 

0.09 

1.586 

3.77 

96 

2x1.875 

231.46 

02/07/2012 

Neptune 

0.05 

0.827 

3.77 

70 

2x1.875 

231.46 

03/07/2012 

Neptune 

0.08 

1.520 

3.77 

79 

2x1.875 

231.46 

04/07/2012 

Titan 

0.10 

1.877 

3.77 

82 

2x1.875 

109.99 

01/08/2012 

Neptune 

0.03 

1.552 

0.4 

53 

4x1.875 


data, the two bands contained 276 baselines, rang¬ 
ing from 21 to 402 m (16 to 309 kA) in band-6, 
and from 21 to 452 m (7 to 167 kA) in band-3. 

In the continuum, we reach a rms of 0.1 
mjy/beam in band-6 and 0.1 mJy/beam in band- 
3. The median rms on the channels associated to 
the , is 6.2 mJy/beam. In the continuum, 

the beam is Oi'Sh x Cf.'60 with a position angle 
(PA) of 88.9° in band-6, and 2"08 x 1"51 with a 
PA of 99.4° in band-3. No emission was detected 
in either the ^^CO(l — 0) or the C^®0(1 — 0) lines 
in band-3. 


2.2. Photometry 


We add two WISE flux at 12 and 22 ^m, two 
Herschel/PACS fluxes at 100 and 160/im, three 
Herschel/Spire fluxes at 250, 350 and 500^m and 
the two mm-fluxes discussed here to the photom¬ 


etry presented by Romero et al. (2012). For the 


Herschel fluxes we adopt the most recent values 
presented by Bustamante et al. (2015). The 


fluxes at wavelengths shorter than 24 /rm were cor¬ 
rected of extinction by applying the dereddening 
relations listed in Cieza et al. (2007). We as¬ 


sume calibration errors of 20% for the optical pho¬ 
tometry and 15% for the photometry up to 24/im. 
The photometry between 100-500 /rm is affected 
by background contamination from the cloud in 
which Sz 91 is embedded, specially at 250, 350 and 
500 /rm (Bustamante et al. 2015). We therefore 


consider the Herschel fluxes at 250, 350 and 500 
as upper limits. For the 870/rm flux we used the 


8.5% uncertainty listed by Romero et al. (2012). 


To derive the fluxes from our ALMA observations 
we perform a fit in the rtri-plane. Because of the 
lack of data at short baselines (i.e., large spatial 
scales), deriving the total flux from the cleaned 
images would result in underestimated values. In¬ 


stead we fit different profiles to the visibilities to 
estimate the flux at the center of the itri-plane. In 
band-6 Sz 91 is resolved (see next section), so we fit 
a Gaussian profile; in band-3 Sz 91 is not resolved, 
so we use a point-source profile. The fluxes com¬ 
prising the SED, their errors and references, are 
listed in Table [U 

3. Continuum and images 

In band-3, Sz91 is only detected (and not re¬ 
solved) in the continuum, with a maximum signal 
to noise ratio (S/N) of « 4.6 and a disk emission of 
0.7 ± 0.1 mjy. The median rms in the individual 
channels is 11.0 mJy/beam. The band-6 results 
are detailed below. 

The continuum image (Fig. left panel) shows 
an inclined disk (the geometry is derived in Sect. 
3.3) with evidences of an inner cavity, as expected 
from the disk’s SED and the previous results of 
T2014. The inner hole is resolved along the pro¬ 
jected major axis of the disk. Integrating over the 
regions of the image with S/N higher than 3 results 
in Fi, 3 m„i = 10.7 ± 1.6 mJy, which is ~ 20% lower 
than the integrated flux estimated from the visi¬ 
bilities (see Table [^. The peak flux in the cleaned 
image is 3.7 ± 0.5 mJy/beam. The difference in 
brightness between the peak flux of the northern 
and southern lobes is below 3tT. The channels cor¬ 
responding to wlsr velocities ranging from 0.5 to 
4.7 km/s show significant (above 3 (t) emis¬ 

sion (Fig. [^. 

Note that the emission at 0.5 km/s, although 
faint when compared to the maximum emission 
of the line, is detected with S/N « 6. The 

corresponding moment 0 and moment 1 images are 
shown in Fig. [^(middle and right panel). The mo¬ 
ment 0 image shows a strong north-south asymme¬ 
try, which is most likely caused by emission from 
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Fig. 1.— Left: band-6 continuum image. Center: Moment-0 image, constructed from the line channels 
showing emission above the Scr level. The white dashed line contours the regions of the image with S/N 
larger than 3. The strong north-south asymmetry is most likely caused by cloud contamination (see text). 
Right: The moment-1 image, constructed using the same channel images as for the moment-0 image. A 
rotating gaseous disk is clearly identified. In all images, north is up, east is left. The white cross shows the 
disk’s center as derived in Sect. 4.3. 


the cloud in which Sz 91 is embedded. Indeed, the 
systemic velocity of the cloud has been measured 
to be ulsr ~ 4.8 km/s (Vilas-Boas et al. 2000 
Tachihara et al. ]|2001 ), which coincides with the 
redshifted emission of Sz91 (southern lobe, 

see the moment-1 image in Fig. right panel). 
The same effect has been already noticed by T2014 
when discussing the (less affected) ^^CO(3—2) line 
profile. These authors estimate that the cloud af¬ 
fects their measurements in the range 4-7 km/s 
■yLSR- Integrating the moment 0 image over the 
regions of the image with emission above 3(T (this 
region is delimited by the white dashed line con¬ 
tour in the central panel of Fig. results in an 
integrated flux of Fi2co = 2.7 ± 0.4 Jy km/s. As 
in the case of the continuum data, this value is an 
underestimation of the true flux due to the clean¬ 
ing process. A Gaussian fit to the visibilities yields 
a larger value of Fi2co = 3.0 ± 0.4 Jy km/s. Be¬ 
cause of the cloud this value still represents a lower 
limit. 

To determine the line center, we used public 
VLT/X-shooter spectra (ID: 089.C-0143) to derive 
the systemic radial velocity of Sz91 (i.e., the cen¬ 
ter of the line) as ulsr ~ 3.4 ± 0.2 km/s, in 
good agreement with the results from Melo (2003) 
and T2014. 


4. Constrains on the star and disk 

In this section, we use our observations in com¬ 
bination with simple modeling to place constrains 
on some of the disk’s properties such as orienta¬ 
tion, geometry, and mm slope. To complete our 
picture of the Sz 91 star-bdisk system we also de¬ 
rive the basic parameters of the central star. 


4.1. Stellar parameters 

First estimations of the Sz 91 stellar parameters 


were obtained by Hughes et al. (1994), who used 
a distance of 140±20 au, an extinction of 2 and 
a spectral type of MO.5. These authors derived a 
stellar mass ranging from M* = 0.49 — 0.69Mq, 
an effective temperature oi T^ff = 3723 K and 
an age ranging from 5-7 Myr. [Romero et 


al. 


(2012) classified Sz91 as a Ml.5 star using high 


resolution spectra, in general agreement with the 
previous result. Using the latter spectral type, 
the most recent distance estimate of ~200 pc, and 
an extinction of = 2, we derive the stellar pa¬ 
rameters using the stellar evolutionary models of 


Siess et al. (2000). We obtain a stellar radius of 


i?* = 1.46i?0, a temperature of Te// = 3720 K, a 
mass of M* = QAIMq, and an age below 1 MyiQ 
In what follows we use these values. 


Siess et al. 


(2000 1 notes that their age estimations for 


stars below 1 Myr tend to be lower when compared with 
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Fig. 2.— Band-6 channels showing significant (above 3cr) emission. The velocity of each channel 
is indicated in the legends. The channel at ulsr = 5.4 km/s does not shows emission, and is shown 

for completeness. The images are plotted in square root strecht. North is up an east is left. The channel 
at ulsr = 3.3 km/s is the closest to the derived systemic velocity of the system (see text). The grey cross 
shows the disk’s center as derived in Sect. 4.3. 


4.2. 


Oi.r] 


slope 


At mm wavelengths, we can use the Rayleigh- 
Jeans approximation to express the flux as oc 
^ Assuming optically thin emission, the mm- 
slope is a function of the dust opacity index /?, 
i.e. Omm = 2 4-/3. In particular, /3 ~ 2 for the 
interstellar medium (ISM) (see [Williams fc Cieza 


2011, and references therein). Trying to fit in this 


way the fluxes at 0.85, 1.3 and 2.7 mm it is not 
possible to match the observed fluxes, as shown 
in Fig. The slope derived from a fit to the 

three fluxes results in oq. 8-2.7mm = 3.36 ± 0.14 
(solid line). Using the two shorter wavelengths 
we obtain ao.g-i.smm = 2.34 ± 0.40 (dotted line), 
whereas using the fluxes at 1.3 and 2.7 mm yields 
to ai. 3 - 2 . 7 mm = 4.07 ± 0.29 (dashed line). This 
large 01 . 3 - 2 . 7 mm would imply moderate or very 
low grain growth. On the other hand, our derived 
ao.8-i.3mm implies /3n.s-i..3mm = 0.34±0.40, which 
suggests the opposite. This latter value agrees 
with the result by T2014 (/3 = 0.5 ± 0.1) and 
with the average values for a set of transition disks 
and protoplanetary disks derived by jPinilla et al. 
[ (2014). In addition, we conclude from our mod¬ 
els (see Sect. 5) that grains of at least 1 mm size 
are needed to match the SED at mm-wavelengths, 
implying that grain growth is happening in Sz 91 
(see Sect. 5). One possible explanation of the 


different slopes is an optical depth effect. En¬ 
hancements in the optical thickness could reduce 
the slope explaining the reduction in ao. 8 - 1 . 2 mm 
compared to the ai. 3 - 2 . 7 mm- However, using a 


sample of 50 disks, Ricci et al. (2012) find that 
low values of a^m can be explained by emission 
from optically thick regions only in the most mas¬ 
sive disks if their sample. Given the relatively low 
mass of Sz91 (Sect. 5), the good agreement with 
the results by T2014, and the significant difference 
between oq. 8 - 1 . 3mm and 01.3-2.7mm, we consider 
that an underestimation of the error bars repre¬ 
sents a more likely explanation of the broken Omm 
slope and that oo.s-i.smm = 2.34 ± 0.40 proba¬ 
bly represents a realistic estimate of the true mm 
slope. 

4.3. Disk orientation 

Given the strong absorption observed towards 
the red-shifted channels, deriving the disk’s center 
and geometry by means of a Keplerian fit to the 


line profile as in e.g. Mathews et al. (20121 would 


estimations from other models. 


result in large uncertainties. Instead, we deter¬ 
mine the disk inclination (j), position angle (PA), 
and center by fitting different disk profiles to the 
continuum visibilities using the Mt/it/MIRIAD 
task. 

A ring morphology matches well our continuum 
observations and the results of T2014. Unfortu¬ 
nately, uvfit doesn’t allow to adjust the thickness 
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Table 2: Photometry of Sz91. All fluxes are 


extinction-corrected. References: (1) Romero et 

al. ( 

2012D, (2) Wright et al. 

(12010 

)• (3) 

Merfn 

et al. 

(2008), (4) Bustamante et al. 

(2015 

1 


Wavelength 

[/rm] 

Flux 

[mJy] 

Error 

[mJy] 

Reference 

0.65 

37.1 

7.4 

1 

1.25 

97.7 

14.7 

1 

1.60 

120.6 

18.1 

1 

2.20 

90.6 

13.6 

1 

3.60 

41.6 

6.2 

1 

4.50 

26.0 

3.9 

1 

5.60 

17.8 

2.7 

1 

8.00 

11.1 

1.7 

1 

12.00 

6.9 

1.0 

2 

22.00 

9.0 

1.4 

2 

24.00 

9.7 

1.5 

3 

70.00 

510.0 

130.0 

4 

100.00 

680.0 

170.0 

4 

160.00 

720.0 

180.0 

4 

250.00“^ 

860.0 

170.0 

4 

350.00’^ 

620.0 

120.0 

4 

500.00“^ 

380.0 

80.0 

4 

850.00 

34.5 

2.9 

1 

1300.00 

12.7 

1.9 

This work 

2700.00 

0.7 

0.1 

This work 


^Considered as an upper limit due to cloud emission. 


of the fitted ring. To test for the effect of using 
different disk profiles, we performed three fits us¬ 
ing a ring, a disk, and a Gaussian profile. Com¬ 
bining these results we derived a PA of 18.2° ± 3.5 
(east of north) and an inclination (i) of 49.5° ±3.5. 
These results are in agreement with those derived 
by T2014 (PA = 17.5° ± 17.7°, i = 40° ± 15°). 
We find the disk’s center to have an offset of 
Aa = -(y.'29 ± O'.'Ol and AS = -O'.'Ob ± 0'.'02 with 
respect to the center of the ALMA field. This off¬ 
set is consistent with the measured proper motion 
of (-20.3,-18.4) mas yr“^ (UCAC3 

2010 |) when including the ^O'.'l astrometric 


al. 


uncertainty from our ALMA Cycle 0 observations 
and the 0706 uncertainty from our input coordi¬ 
nates (2MASS catalogue, Cutri et al. [2003 ). We 
corrected for this offset in our analysis of the disk. 



1000 


Wavelength [/im] 


Fig. 3.— Different fits to the mm-slope (omm)- It 
is not possible to simultaneously match the three 
measurements (see text). 


4.4. Outer disk 

A quick inspection to the continuum (Fig. 
left panel) and the ^^CO moment-0 and moment- 
1 (Fig. [1 central panel) images shows that the 
gaseous disk extends further out than the contin¬ 
uum disk. The resolution of our images along the 
major axis of the disk is set by the beam size in 
that direction i.e. 0761. Keeping this in mind, 
we produced a 3-pixel wide cut along the blue- 
shifted (i.e., less affected by cloud emission) semi¬ 
major axis of the disk for the continuum and ^^CO 
moment-0 images (Fig. |^. We compute the fluxes 
every Cf.'16, which means that they are correlated 
(i.e., their separation is smaller than the beam 
size). Still, we can use them to obtain informa¬ 
tion about the extension of the gaseous and dusty 
disks. In the plot, the vertical dotted lines indicate 
the position at which the emission becomes signif¬ 
icant (i.e., above 3cr). The continuum emission 
becomes undetectable beyond ~ 171 (i.e., ^220 
au) from the center, whereas the ^^CO emission 
is still detected at ~ 2” (^400 au), i.e., almost 
1.5 beams further out than the dusty disk. We 
therefore conclude that, although we are limited 
by the moderate resolution of the beam size, the 
gaseous disk is undoubtedly detected at larger dis¬ 
tance from the star than the dusty disk. The outer 
radius of the dusty and gaseous disk estimated by 
T2014 at 345 GHz (170 ±20 au for the continuum, 
420 au for the gas) matches notably well with our 
results. Our limited spatial resolution doesn’t al¬ 
low us to probe for the sharp outer edge in the 
continuum disk observed in other disks, e.g., TW 
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Fig. 4.— Cut along the blue-shifted (i.e., not af¬ 
fected by cloud extinction) semi-major axis of the 
disk. The points are correlated (i.e., their sep¬ 
aration is smaller than the beam size along the 
disk’s major axis). The vertical black and red dot¬ 
ted lines show the position at which the emission 
is higher than 3x the rms of the continuum and 
moment-0 images, respectively. The nega¬ 

tive values are likely artifacts from the cleaning 
algorithm. 


the presence of a cavity (see Hughes et al. |[20^ 
Andrews et al. 1[MI] iHughes et al. ||2009D. Witt 


Andrews et al. 1[MI] IHughes et al. ||2009D . With 

our observations we cannot estimate the sharpness 
of the inner wall (see e.g., Mulders [2013 ). In the 
next section we derive the cavity size in the con¬ 
tinuum by fitting the uv plane. 

The visibilities have lower S/N than the 

continuum and they are severely affected by the 
cloud. This complicates the description of the 
gaseous disk in terms of the analysis of the vis¬ 
ibility profile. However, if the gaseous disk veloc¬ 
ities are governed by Keplerian rotation, we can 
use the line profile to estimate how close to 
the central star we detect gas using the projected 
observed velocity v(r)Kep = sin b As¬ 

suming that the line center is at ?;lsr = 3.4 km/s, 
we detect emission down to ulsr = 0.5 km/s 
with a S/N of ^6 (Fig. [^. This corresponds to 
—2.9 km/s in the star’s reference system. Using 
this velocity, the stellar mass and inclination pre¬ 
viously derived, and a distance of 200 pc, we con¬ 
clude that our observations probe the gas 

down to ~ 28 au from the central star. 


Hya (Andrews et al. 20111 and HD 163296 (de 


Gregorio-Monsalvo et al. 


20131 that could be a 


consequence of radial drift (Birnstiel et al. |2013D . 


4.5. Inner disk 

The shape of the continuum cleaned image is 
indicative of the presence of an inner cavity. To 
explore the partially resolved inner disk we here 
focus on the radially averaged deprojected ALMA 
visibilities. We split the continuum and vis¬ 

ibilities into 30 radial bins, ensuring that all bins 
contain the same amount of visibilities. This re¬ 
sults in minimum and maximum bin sizes of 3.7 
kA and 22 kA, respectively, and represents a good 
compromise between uv sampling and error bars. 
The real part of the binned visibilities of our band- 
6 observations as a function of uv distance are 
shown in Fig. The steeper slope of the 
visibilities shows that the detected gaseous disk is 
indeed more extended than the dusty disk, as al¬ 
ready noticed from the cleaned images. The null in 
the continuum visibilities at ^ 150kA indicates a 
discontinuity caused by the lack of mm-sized par¬ 
ticles in the inner disk regions, or, in other words. 


5. Radiative transfer model 


Based on the obtained observational con¬ 
straints, in what follows we develop an axisym- 
metric model for the disk around Sz 91 using the 


radiative transfer code MCFOST (Pinte et al. 


2006, 20091. MCFOST computes the dust tem¬ 


perature structure under the assumption of radia¬ 
tive equilibrium between the dust and the local 
radiation field. Provided the proper set of param¬ 
eters, MCFOST creates synthetic SEDs as well as 
continuum and ^^CO images. From these images 
we create synthetic interferometric observations 
with the same uv coverage as our observations, 
which are compared with our observations. 

We start with a description of our model as¬ 
sumptions about the dust structure, composition, 
and distribution. Then we use a genetic algorithm 


(Charbonneau 1995) to identify a best-fit model 


that can reproduce our observations. 


5.1. Dusty Disk Structure 

To account for the SED shape at mid-infrared 
wavelengths we define two inner disk-regions: the 
innermost region extending from the dust subli¬ 
mation radius to Rout,i which is severely dust- 
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Fig. 5.— Deprojected, radially binned real and imaginary visibilities as a function of uv distance for the 
continuum (left) and the gas (right). For the continuum, each bin shows the sum of the visibilities averaged 
over the continuum channels. For the , each bin shows the sum of the visibilities over the channels 

containing signihcant emission (see text). The null in the continuum visibilities at ^ 125kA indicates the 
presence of a large inner cavity in the mm-sized dust distribution. The red solid lines show the model 
discussed in Sect. 5. 


depleted (“Region 1”) and accounts for the fluxes 
at 12 /rm and below, and the intermediate region 
extending from Rin ,2 to Rout ,2 (“Region 2”) where 
some dust is emitting at 22 and 24 /rm. There is 
no observational constrain on the surface density 
distribution of these two regions, and therefore we 
adopt a simple power-law of the form 


where Sioo is the surface density distribution 
at 100 au and Ps is the power-law index. We de¬ 
fine the outer disk, “Region 3”, as the part of the 
disk extending from the cavity radius (i?cav) out¬ 
wards. For a number of disks where the gaseous 
disk is detected beyond the dusty disk as in Sz91 
a tapered-edge profile can naturally explain the 


observations (e.g.. 

Hughes et al. 

to 

o 

o 

00 

Williams 

ik Cieza 201 1[ Andrews et al. | 

2011 

Mathews 

et al. 

2012| [Cieza et al. 2012| Zhang et al. 

20141. 

There are three exceptions to this behav- 


ior. To explain the high resolution observations 
of IM Lup, TW Hya and V4046 Sgr strong vari¬ 
ations of the gas to dust ratio at large radii are 


required 1 

Panic et al. 

2009 Andrews et al. |2012 

Rosenfeld et al. |2013 

. Our data do not have the 


high spatial resolution and sensitivity needed to 
test for radial changes in the gas to dust ratio and 
therefore we adopt the tapered-edge prescription 
to describe the outer disk. The surface density 
distribution in Region 3 is then defined as 


E(r) = Sc f ^ exp 



( 2 ) 


where Rq is the characteristic radius which de¬ 
fines where the tapered-edge begins to dominate 
over the power-law component. Sc is the surface 
density at Rc, and 7 corresponds to the viscosity 
power-law index in accretion disk theory (zz oc R'^, 
Lynden-Bell fc Pringle ||1976[ [Hartmann et al. | 
19981. For each region we define the scale height 


H{r) assuming that the dust follows a Gaussian 
vertical density profile 


H{r) = i?o(r/100au)^. 


(3) 


where tjj is the flaring parameter of the disk and 
Hq is the scale height at r = 100 au. 

5.2. Dust composition and distribution 

We assume homogeneous spherical dust parti¬ 
cles. The scattering opacities and phase functions, 
extinctions and Mueller matrices are computed 
using Mie theory. For the composition we as¬ 


sumed compact astro-silicates (Mathis & Whiffen 


1989 1 with a power-law size distribution (dn{a) oc 
a~Pda, with p = 3.5 e.g., Draine [[2006 ). This as¬ 
sumption requires minimum and maximum grain 
sizes of ttinin = 2 /rm,aniax = 15/xm in region 1 to 
not overproduce the flux at 8 and 12 pm due to 
the silicate resonance bands. We note that using 
a different grain composition as e.g. carbon grains 
it is also possible to match the SED using smaller 
values of Omin inside the cavity. In any case, the 
lack of emission above photospheric levels below 
8 /im indicates that the innermost region of the 
cavity must be severely depleted of dust particles. 
Region 2 is mainly constrained by the fluxes at 




























































































22, 24 and 70 fj,ia and therefore many combina¬ 
tions of aniin, amax, and the dust mass can repro¬ 
duce the observed fluxes. The minimum grain size, 
however, cannot exceed Omin = 0.1/im in order to 
match the observed fluxes. 


5.3. Fitting procedure 


We use the stellar parameters derived in Sect. 
4.1 and the inclination {{) and position angle (PA) 
derived in Sect. 4.3. We focus on identifying a 
model that can describe the thermal structure of 
the disk (constrained by SED and continuum ob¬ 
servations) using the disk structure and dust pa¬ 
rameters described in the previous sections. Re¬ 
gions 1 and 2 are poorly constrained and therefore 
we fix the power-law indexes of their surface den¬ 
sity distribution to the standard value of = 1 


following Andrews & Williams (2007). In Region 


1 we fix the minimum and maximum grain sizes to 
the values discussed in the previous section. Using 
this prescription we find that we need a maximum 
dust mass of Mdust,i = 1 x 10“^ to reproduce 
the SED shape below 12/rm. We adopt this value 
to describe the dust content inside Region 1. A 
first exploration of the Hq parameter in this re¬ 
gion reveals that it is completely unconstrained, 
and we adopt an arbitrary value of ifoi = 5au. 
In Regions 1 and 2 we fix the minimum grain size 
to amin = 0.05/rm, as in the interstellar medium. 
For Region 3 we explore the effects of varying the 
surface density parameter 7 within the range of 
[—1,1]. We find that our data is equally well fit¬ 
ted with values ranging from 7 = [—0.2, 0.8] and 
use 7 = 0.3. Similarly, our data does not pro¬ 
vide observational constraints on the flaring index 
(i/)) in any of the three regions. For all three re¬ 
gions we adopt a typical value of ip = 1.15 from 
model results in , e.g. T Cha (Cieza et al. [20111, 


IM Lup (]Pinte et al. 2008) and HD 163296 (de 
Gregorio-Monsalvo et al. 120131 ). Free parameters 


in our model approach are the outer radius of Re¬ 
gion 1 and 2 {Routi^ 2 )^ the inner radius of Region 2 
(Rm,2), and the cavity and characteristic radius of 
Region 3 {Rcav,Rc)- We also fit the dust masses 
Mdust^ 3, the maximum grain sizes Umax^ 3 and the 
scale heights i7o2,3- 

We explore the parameter space by means of a 
genetic algorithm. The procedure starts by ran¬ 
domly selecting a number of models. The param¬ 
eters of these models are defined as their “genes”. 


Models producing a better match to our obser¬ 
vational dataset have the highest chance to re¬ 
produce their parameter values (genes) into the 
next generation of models. The quality of a fit is 
measured with the combined reduced given by 
the sum of reduced from SED and the contin¬ 
uum visibilities, i.e = Xsed + xlv We use 
a population of 50 “parent” models and let them 
evolve through 50 generations, resulting in a total 
of 2500 models. We find that after 35 generations 
the models quickly converge, and after 45 gener¬ 
ations the variations of each free parameter are 
within the 5% of the average value of the param¬ 
eter. 

5.4. Modeling results 

We find a best-fit model that agrees well with 
our ALMA observations and with the SED of 
Sz91. Its parameters are summarized in Table 
The most significant result of our modeling is that 
a large cavity of ^ 97 au in the dust grain distri¬ 
bution is required to mach our observations. This 
is clearly illustrated by the null in the real part of 
the continuum deprojected visibilities as shown in 
Fig.(left panel). For direct comparison with our 
ALMA band-6 images we created synthetic ALMA 
observations using the SIMALMA/cas.^ package. 
We created an antenna configuration file that 
mimics our observations using the buildConfig- 
urationFile/Analysis Utilitie^ task. We 
then run the simulation using the same PWV, ex¬ 
posure time and hour angle as in our observations. 
The synthetic ALMA image of our model for Sz 91 
and the residuals obtained by subtracting the syn¬ 
thetic image from the real observations are shown 
in Fig. (central and right panels). The residual 
image has an rms of 0.1 mJy/beam similar to the 
rms of the observations. Interestingly, our simu¬ 
lated image also shows a brightness difference be¬ 
tween the northern and southern lobes of the disk, 
with the northern lobe being the brighter one (al¬ 
though our MCFOST model image is axisymmet- 
ric). To investigate this feature we repeated the 
SIMALMA/CASA simulations using different uv 
coverages and PWV. We find that the asymmetry 


®More information about the CASA package can be found 
at: http://casa.nrao.edu/ 

^Analysis Utilities is an external package provided by NRAO 
to complement the CASA utilities. It can be found at: 
http: //casaguides.nrao.edu/index.php?title=Analysis_Utilities 
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Fig. 6.— Left: Band-6 continuum observations. Center: Synthetic, ALMA simulated images. Right: 
Residual image. Notice the different color scale used in this figure. 


is a consequence of the low S/N of the visibilities 
at longer baselines as increasing the uv coverage 
at larger baselines or decreasing the PWV of the 
simulated ALMA observations removes the appar¬ 
ent difference in brightness. 

The SED and surface density distribution of the 
model are shown in Figs.[^and|^ respectively. In 
Fig.[7] we additionally show the individual contri¬ 
butions of Region 1 and Region 2. 

5.5. Model uncertainties 

To estimate which model parameters are well 
constrained from our observational dataset we pro¬ 
ceed as follow. We fix all but one parameter and 
run a grid of models exploring the local param¬ 
eter space. We then derive the reduced of 
the models forming this local grid. This proce¬ 
dure is repeated for all free parameters. Keep¬ 
ing in mind that this is equivalent to assume that 
all parameters are uncorrelated (which is clearly 
not the case), this approach allow us to get a first 
idea of which parameters are clearly constrained 
by our observations without spending unreason¬ 
able amounts of computing time. The results of 
this exploration are shown in Fig. These plots 
illustrate the severe degeneracy of most of the pa¬ 
rameters describing Regions 1 and 2 in our model 
disk. In contrast, we find that the parameters such 
as Mdust, amax and the cavity size {Rcav) in Re¬ 
gion 3 are more constrained by our observations. 
In particular, values of Rcav below 90 au are very 
difficult to reconcile with the ALMA observations. 
According to our model, most of the disk’s mass 



Fig. 7.— Observed SED of Sz91, in red symbols. 
The simulated star’s photosphere is shown in grey. 
The model described in Sect. 5 is shown as a solid 
black line. The dotted and dashed lines represent 
the separate contributions of Region 1 and 2 of 
our model disk, respectively. 

(10 M0,~ 70% of Mdust, 3 ) is concentrated in a 
ring ranging from 97 to 140 au. 

5.6. Gaseous Disk 

Our observations trace the optically thick ^^CO 
line which is severely affected by the cloud. We 
therefore only present a simple model of the gas 
component of the disk around Sz 91 that is able 


10 
























E [gr/cm^] 



Fig. 9.— (i-e., the difference between the reduced ^ind its minimum) of the fixed parameters in 

our modeling (see text). Blue dashed line corresponds to Region 1, red-dotted and black-solid lines indicate 
Regions 2 and 3, respectively. 



Fig. 8.— Surface density distribution of the 
model. The dotted and dashed lines represent 
Rin ,2 and Rcav, respectively. The model extends 
down to 0.025 au (sublimation radius). 


ing the gas to dust ratio and total gas mass as in 


Williams & Best (2014). We fix the gas to dust ra¬ 


tio in the outer region (Region 3) to the standard 
value of 100, and we ensure that the freezes- 
out at 20 K ( Qietal. |2004[ de Gregorio-Monsalvo 
et al. We assume that the gas and the dust 

have the same temperature (local thermal equilib¬ 
rium, LTE). For the abundance we use the 

standard H2 / ratio (10"^), which is a reason¬ 

able value for protoplanetary disks as confirmed by 


France et al. (2014|. The gas is assumed to be in 


Keplerian rotation with an inner radius of at least 
28 au (in agreement with our observations). As¬ 
suming that the gas is in hydrostatic equilibrium 
in this direction, the gas profile is described by 


p(r, z) = p(r, 0)exp 


2H{ry 


(4) 


to match the line profile reasonably well us¬ 

ing standard assumptions. We only consider the 
blue-shifted (less absorbed) side of the line. Given 
the sparse observational constrains on the gaseous 
disk, we build this model upon the best-ht model 
for the dusty disk obtained in the previous sec¬ 
tion. The lack of measurements of the isotopo- 
logues ^^GO and G^®0 prevents us from estimat- 


The scale height of the gas component H{r) is 
given by 


H{r) 


I r^kBT{r) 


(5) 


where fee is the Boltzmann constant, G is the grav¬ 
itational constant and mn is the Hydrogen mass. 
We also take into account the thermal broadening 
of the velocities (defined as Vth = \/2/cbTco/wcO) 
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Table 3: Model Results. Fixed parameters are 
listed above the line._ 


Parameter 

Value 

T* 

3720 K 

R* 

1.46 Rq 

M* 

0.47 Mq 

d 

200 pc 

i 

49.5° 

PA 

18.2° 

P 

3.5 

^mini 5 ^maxi 

2,15/im 

^min2 3 

0.05/rm 

Ho, 

5 au 

Ps,a 

1 

7 

0.3 

Vl,2.3 

1.15 

Afdusti 

1 X 10-4 M® 

Ri 

0.025 au 

^max2 

5/im 

^max3 

1000/rm 

Rout\ 

85 au 

^in2 

85 au 

^out2 

97 au 

Rcav 

97 au 

Rc 

100 au 

'^^dust2 

0.7 M® 

-^dusts 

14.3 M® 

Ho, 

5 au 

Ho, 

10 au 

log(Siooi) 

—8.20 [gr/cm^]’^ 

log(Sioo2) 

—2.95 [gr/cm'^]^ 

log(Uc) 

-1.98 [gr/cm^]’^ 


“Derived values. 


'^C0(2-1) 



Velocity (km/sec) 

Fig. 10.— Observed line profile, in black his¬ 
togram. The velocity scale is centered at the ref¬ 
erence system of the star, using ulsr = 3.4 km/s. 
The profile was computed by integrating the flux 
contained in the > 3 ct region of the moment-0 
image (white dashed contour in central panel of 
Fig.§ in each channel. The vertical dotted lines 
show the center of the line and the highest velocity 
bin observed at > 3 ct significance, corresponding 
to -2.9 km/s. The vertical dotted line at 2.9 km/s 
highlights the strong absorption-like effect on the 
red-shifted side of the line. The red his¬ 

togram shows the gas model presented in Sect. 5. 

and the surface density of this model are shown 
in Fig. (right panel) and Fig. (red line), re¬ 
spectively. The integrated line profile of the 
model, along with the real observations, is shown 
in Fig. 

6. Discussion 


where mco is the molecular mass of the ) 

and the effect of turbulence (r^turb) on the 
emission along the line of sight. 

Under these assumptions, we manually explore 
different models by varying the turbulence veloc¬ 
ity (t'turb) within a range of [10,300] m/s (in accor¬ 
dance with the observations discussed by [Hughes] 
) and the gas mass inside the dust 
depleted regions (Region 1 and 2). We find a 
model that agrees relatively well with our obser¬ 
vations using a continuous gaseous disk (i.e. with¬ 
out depletion or inner cavity) with Uturb = 120 
m/s. The radially binned deprojected visibilities 


6.1. Are three zones required? 

Assuming three different disk regions we find a 
model that reproduces reasonably well our obser¬ 
vational dataset. As a separate exercise we also 
tried to fit a disk model composed by only two 
regions (the inner cavity and the outer disk). Af¬ 
ter running the genetic search we could not find a 
good match to our observations with such a disk 
structure. This result is in agreement with the 
findings of T2014, who need an extra component, 
either a circumplanetary disk or a narrow, opti¬ 
cally thick inner ring, to reproduce the fluxes at 
22 and 24/im. 
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6.2. Model limitations 

The local parameter exploration (Fig. repre¬ 
sents a useful approach to estimate how sensitive 
are the model parameters to our observations. We 
find that most of the parameters describing Region 
1 and 2 remain unconstrained from our observa¬ 
tions. For instance, in our representative model 
there are no empty gaps as we find Rout,i = Rin ,2 
and Rout .2 = Rcav, although other values for these 
parameters (i.e., an empty gap between Region 1 
and 2 , or a dust empty innermost region) could 
equally reproduce our observations (as indicated 
by the profile in the Rin and Rout panels in 
Figj^. For all 3 regions we find that the the expo¬ 
nents of the surface density distributions are com¬ 
pletely unconstrained by our observations. The 
same happens with the flaring index tjj, as ex¬ 
pected if the emission is optically thin at these 
wavelengths. However, the exercise we performed 
does not allow to derive proper statistical errors 
of the fitted parameters as we did not investigate 
correlations between the parameters. The latter 
would require to map the entire multidimensional 
parameter space which is beyond the scope of this 
paper and which would require to spend an un¬ 
reasonable amount of computing time. We re¬ 
serve this exercise for the future, when data with 
more sensitivity and much better angular resolu¬ 
tion will be available. Regarding the gas model, 
we find that it is difficult to reproduce the ob¬ 
served profile with our model. Overall we 

see that compared to the observations the model 
shows a lack of emission at higher velocities, and 
an excess at lower velocities. Even using a non- 
depleted gaseous disk our model produces less flux 
at higher velocities than the observations. We note 
that this does not exclude the possibility of a gas- 


depleted cavity (as found for HD 142527, Perez et 


al. 20141, since we also find that it is possible 


to produce a similar profile modeling a gas- 

depleted cavity using higher values of turbulence, 
or using different surface density profiles for the 
gaseous and dusty disks (i.e., different 7 for the 
gas and dust disks). 

6.3. Cavity-making mechanism 

With a size of 97 au in radius, Sz 91 has by far 
the largest the inner cavity observed in a transi¬ 
tion disk around a T Tauri star. This cavity size 


and accretion rate (7.9 x 10 “^^Mq yr“^) exclude 
photoevaporation alone as the mechanism respon¬ 
sible for the cavity (e.g., Rosotti et al. | 2013 


Owen et al. 2012). Similarly, studies of grain 


growth conclude that grain growth alone cannot 
produce the large cavities resolved by interfero¬ 


metric mm observations (Birnstiel et al. 2013). 


Romero et al. (2012) excluded a stellar compan¬ 
ion down to separations of ~ 30 au with a contrast 
of AK ~ 3.3. Using NextGen models (e.g. Allard 


et al. 1997), an age of 1 Myr and the stellar 


mass derived in Sect. 4.1 translates into a sensi¬ 
tivity close to the brown dwarf limit of ^ 80 Mj. 
In addition, Melo (2003) found no evidence for 


a close-in binary companion around Sz 91 in their 
3-years radial velocity survey (down to masses of 
« O. 2 M 0 ). However, we cannot yet exclude that 
a low mass stellar companion is carving the cav¬ 
ity in the disk around Sz91. The other possible 
mechanism is planet formation. In fact, our direct 
detection of gas inside of the large dust-depleted 
cavity, the presence of /rm-sized dust inside the 
cavity, and the large cavity size indicated by the 
continuum visibilities matches well several predic¬ 
tions of models for multiple (Jupiter-sized) planet 
formation (e.g. Lubow fc D’Angelo |2006[ Dodson- 


Robinson fc Salyk ||2011 Zhu et al. ||20li ). 


6.4. Inner disk structure: planet forma¬ 
tion signature? 

There are at least four transition objects with 
detectable (but different) accretion rates that can 
be explained by an inner disk structure simi¬ 
lar to the one derived for Sz91. DMTau and 
RXJ1615.3-3255 have resolved cavities and re¬ 
quire an empty disk very close to the star be¬ 
cause of the lack of NIR excess, but some dust in 
the outer regions of the cavity to account for the 


MIR emission (Andrews et al. 2011). However, 


a direct comparison with Sz 91 is difficult because 
the modeling approach adopted by [Andrews et alT 
I ( [2011 ) is different from ours. RX J1633.9-2442 
and [PZ99] J160421.7-213028 have been analyzed 
using similar tools than those used in this work 
for Sz91. In both systems, an empty inner disk 


and some dust inside the main cavity (Cieza et 


and 

some 

al. 

2012 

2014 

) are 


2012 Mathews et al. 2012[ [Zhang et al. 


2014) are needed to explain the mm-observations 


and SED. Interestingly, the mm-sized dust distri¬ 
bution in [PZ99] J160421.7-213028 is concentrated 
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in a ring similar to what we find for Sz91 (Zhang 


et al. 2014 Mathews et al. 2012). Pinilla et 


al. (2012) concludes that this ring-shaped dusty 


structure is a product of the planet-disk interac¬ 
tion. 

While in all five cases strong parameter degen¬ 
eracies of the models do not allow us to constrain 
the exact configuration of the dust inside the cav¬ 
ity, photoevaporation and grain growth can be 
ruled out as the mechanisms carving the inner 
hole of these disks. Furthermore, for DMTau, 
RXJ1633.9-2442, and [PZ99] close binarity can 
also be excluded and planet-formation has been 
suggested to explain both the cavities and the 


presence of dust inside the cavity (e.g.. 

Rice et al. 

2006 

Cieza et al. 

2012 

Mathews et al. 

2012 ). 

As suggested by 

Cieza et al. 

(2012), these inner 


structures are likely simple approximations to the 
complex structures predicted by hydrodynamical 


models of planet formation (Dodson-Robinson & 


Salyk 2011 

, which have now been directly ob- 

served (e.g.. 

Casassus et al. 

2013 Quanz et al. 


2013) in a few systems. 


7. Conclusions 


We presented ALMA band-6 and band-3 obser¬ 
vations of Sz 91 which clearly show that there must 
be a very large inner cavity in the disk around 
Sz91. We also detect gas inside the cavity 

down to at least 28 au. We detect gas at larger 
radii (400 au) than the dust (220 au). We used a 
genetic algorithm to construct an MCFOST model 
of the disk around Sz 91 explaining our ALMA ob¬ 
servations and the SED of the system. We find 
that a three-zone model is required to fully ex¬ 
plain the observations. The inner region of the 
disk is significantly depleted of dust particles, and 
there are « O.7M0 of dust confined within the 85- 
97 au region. The exact spatial distribution of the 
dust in these two regions is not constrained by the 
available observations but our uu-fit to the data 
indicates that the dust depleted region extends to 
~ 97 au. This implies that there must be signifi¬ 
cant amounts of gas inside the dust depleted inner 
regions. The bulk of the dust mass is located in 
a ring extending from 97 au to 140 au. The dif¬ 
ferent outer radius observed for the and the 

dust can be explained by an exponential tapered 
surface density profile. 


A yet undetected, very low mass companion 
(below O. 2 M 0 ) in the inner 30 au, could create the 
large cavity observed in Sz 91. However, the size of 
the disk’s cavity and the presence of gas and small 
dust particles inside the cavity agree remarkable 
well with the theoretical predictions for multiple, 
Jupiter-sized, planet formation. Considering its 
spectral type and low stellar mass, Sz91 has an 
extremely large inner cavity when compared with 


other transition disks (e.g. Andrews & Williams 


2007). Finally, we find that similar structures of 


the inner regions of dusty disks are also observed in 
at least four other planet-forming candidates tran¬ 
sition disks. It is tempting to interpret the inner 
disk structure and large cavity size as a signpost 
of forming planets inside these disks. 
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